/*******************************************************************************

This code file produces Panel B of Figure A24, "421-a Participation Under Alternative Tax Incentives and Inclusionary Shares."

*******************************************************************************/

*** Manage settings

	run "$dir/code/modules/settings.do"

********************************************************************************
* "Laffer curve" analysis of inclusionary housing
********************************************************************************

	set seed 1234

*** Bootstrap estimation of Laffer curve

	forvalues b = 1(1)1001 {
		
		use "$data/clean/cleaned_data.dta", clear
		bsample, cluster(nta)
		
		quietly summ dtaxrate_onsite
		local avgrate = r(mean)
		
		capture drop xb_pred 
		quietly logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		predict xb_pred, xb
		
		mat b = e(b)
		local coef_dtaxrate = b[1,1]

		forvalues share = 0(0.005)0.5 {
		
			capture drop pr_pred sh_units_inclusionary
		
			local int_share = int(round(100*`share',1))
			
			preserve
						
			quietly gen pr_pred = exp(xb_pred-`coef_dtaxrate'*dtaxrate_onsite + `avgrate'*`coef_dtaxrate'/`share')/(1+exp(xb_pred-`coef_dtaxrate'*dtaxrate_onsite + `avgrate'*`coef_dtaxrate'/`share'))
			
			quietly egen inclusionaryunits = sum(pr_pred * floor(`share' * unitsres))
			quietly egen allunits = sum(unitsres)
			
			quietly summ inclusionaryunits
			local numerator = r(sum)
			
			quietly summ allunits
			local denominator = r(sum)
			
			local outcome = `numerator' / `denominator'
	
			capture use "$data/clean/boot_laffer.dta", clear
			local ob = 100*(`b'-1) +  `int_share'
			quietly replace b = `b' if _n == `ob'
			quietly replace sh_inclusionary = `share' if _n == `ob'
			quietly replace sh_units_inclusionary = `outcome' if _n == `ob'
			
			di "Observation: `ob'"
			di "Bootstrap: `b'"
			di "Requirement: `int_share'"
			di "Outcome: `outcome'"
			
			quietly save "$data/clean/boot_laffer.dta", replace
			
			restore
		
		}
		
	}
	
*** Save ground truth

	use "$data/clean/cleaned_data.dta", clear
	
	summ inclusionary_onsite [aw=unitsres]
	
	local sh_inclusionary = 0.2
	local sh_units_inclusionary = 0.2 * r(mean)
	
*** Chernozhukov et al procedure for simultaneous confidence intervals

	use "$data/clean/boot_laffer.dta", clear
	
	replace sh_units_inclusionary = 0 if sh_inclusionary == 0
	
	bys sh_i: egen sh_units_mean = mean(sh_units_inclusionary)
	bys sh_inclusionary: egen sd = sd(sh_units_inclusionary)
	gen t = abs(sh_units_inclusionary - sh_units_mean)/sd
	bys b: egen max_t = max(t)
	summ max_t, d
	
	local p95 = r(p95)
	di `p95'
	
	gen lo = sh_units_mean - `p95'*sd
	gen hi = sh_units_mean + `p95'*sd
	
	sort sh_inclusionary
	
*** Produce graph

	tw (line sh_units_mean sh_inclusionary if b==2, lcolor(edkblue)) || ///
		(rline hi lo sh_inclusionary if b==2, lcolor(ltblue)), ///
		xline(`sh_inclusionary', lpattern(dash) lcolor(gs9)) ///
		yline(`sh_units_inclusionary', lpattern(dash) lcolor(gs9)) ///
		graphregion(color(white)) ylabel(, nogrid) ///
		legend(order(1 "Bootstrap Mean" 2 "95% Simultaneous Confidence Band")) ///
		ytitle("Citywide Inclusionary Share of All Units") xtitle("Set-Aside Share")
		
	graph export "$figs/citywide_laffer_curve_nofe.pdf", replace
	graph export "$figs_overleaf/citywide_laffer_curve_nofe.pdf", replace
	
/*
	
*** Find argmax (inclusionary share peak of Laffer curve)

	forvalues b = 1(1)1001 {
		
		use "$data/clean/cleaned_data.dta", clear
		bsample, cluster(nta)
		
		capture drop xb_pred 
		quietly logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		predict xb_pred, xb
		
		mat b = e(b)
		local coef_dtaxrate = b[1,1]

		forvalues share = 0(0.001)0.20 {
		
			capture drop pr_pred sh_units_inclusionary
		
			local int_share = int(round(100*`share',0))
			
			preserve
						
			quietly gen pr_pred = exp(xb_pred-`coef_dtaxrate'*dtaxrate_onsite + `avgrate'*`coef_dtaxrate'/`share')/(1+exp(xb_pred-`coef_dtaxrate'*dtaxrate_onsite + `avgrate'*`coef_dtaxrate'/`share'))
			quietly gen sh_units_inclusionary = pr_pred*`share'

			quietly summ sh_units_inclusionary [aw=unitsres]
			local outcome = r(mean)
	
			capture use "$data/clean/boot_laffer_2.dta", clear
			local ob = 100*(`b'-1) +  `int_share'
			quietly replace b = `b' if _n == `ob'
			quietly replace sh_inclusionary = `share' if _n == `ob'
			quietly replace sh_units_inclusionary = `outcome' if _n == `ob'
			
			di "Observation: `ob'"
			di "Bootstrap: `b'"
			di "Requirement: `int_share'"
			di "Outcome: `outcome'"
			
			quietly save "$data/clean/boot_laffer_2.dta", replace
			
			restore
		
		}
		
	}


	bys b: egen max_sh_units = max(sh_units_inclusionary)
	bys b: egen max_sh_inclusionary = max(sh_inclusionary*(sh_units_inclusionary == max_sh_units))
	
	keep b max_sh_inclusionary
	duplicates drop

*/
